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Abstract. Full relativistic simulations in three dimensions invariably develop 
runaway modes that grow exponentially and are accompanied by violations of the 
Hamiltonian and momentum constraints. Recently, we introduced a numerical method 
(Hamiltonian relaxation) that greatly reduces the Hamiltonian constraint violation 
and helps improve the quality of the numerical model. We present here a method 
that controls the violation of the momentum constraint. The method is based on the 
addition of a longitudinal component to the traceless extrinsic curvature , generated 
by a vector potential to,, as outlined by York. The components of Wi are relaxed to 
solve approximately the momentum constraint equations, pushing slowly the evolution 
toward the space of solutions of the constraint equations. We test this method with 
simulations of binary neutron stars in circular orbits and show that effectively controls 
the growth of the aforementioned violations. We also show that a full numerical 
enforcement of the constraints, as opposed to the gentle correction of the momentum 
relaxation scheme, results in the development of instabilities that stop the runs shortly. 
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1. Introduction 

For most of the past decade, the main theoretical thrust in gravitational research has 
been directed toward obtaining stable and accurate numerical models of compact-object 
binary systems. One of the most difficult problems to tackle has been the control 
of exponentially growing instabilities that degrade the quality of any simulation and, 
eventually, terminate it. General relativistic numerical simulations of BNS systems have 
made significant progress in the past years [U El EH IH 03 HE E| , in particular with the 
addition of advanced numerical techniques such as AMR [S]. 

The numerical formulations employed to simulate astrophysical systems can be 
divided in unconstrained and constrained methods. Unconstrained formulations like 
ADM and BSSN [TQl IH] simply evolve the gravitational fields without any attempts 
at controlling the violation of the time independent Hamiltonian and momentum 
constraints; these equations are merely monitored to gauge the accuracy of the model. 
Constrained formalisms, on the other hand, enforce the satisfaction of the constraint 
equations at either the analytical level (i.e., ingraining them in some way in the time 
evolution equations) or at the numerical level (i.e., regularly solving the constraints 
during the evolution). While constrained methods have been applied mostly to problems 
with spherical 12 and axial symmetry ^3JEJE3, there have been recent applications 
to three-dimensional scenarios [TH} I17j. 

We propose in this and a companion article JH] (from now on, Paper I) an 
alternative type of evolution in which the Hamiltonian and momentum constraints 
are only approximately solved at every time step, gently steering the evolution toward 
the space of solutions of these equations without completely forcing their numerical 
satisfaction. In Paper I we described a method that controls the Hamiltonian constraint 
violation (Hamiltonian relaxation or HR). In this paper we introduce a complementary 
scheme that reduces the momentum constraint residuals (momentum relaxation or MR). 
The constraint relaxation techniques utilize the conformal decomposition of the spatial 
metric and extrinsic curvature presented by York [HIl ED] , which has traditionally been 
used to solve the initial value problem for binary systems [2H 122 EH1 EH 125] • In 
this decomposition, a conformal factor ip factored out of the spatial metric and a 
longitudinal addition to the extrinsic curvature generated from a vector potential Wi are 
used to satisfy the Hamiltonian and the three components of the momentum constraint 
respectively. HR drives ip to the solution space of the Hamiltonian constraint by means 
of a parabolic equation for the conformal factor, in the spirit of the K-Driver [2^] used for 
the lapse. The MR method described in this paper uses Wi to push the simulation toward 
the space of solutions of the momentum constraint. In both full relaxation of tp 

and Wi would lead to the numerical solution of the constraints. However, the stability 
of the relaxation methods relies on gently updating these fields during the evolution. 
We show in Appendix A that a full relaxation scheme becomes unstable rather quickly 



when used in combination with BSSN. 

We tested the new algorithms by simulating binary neutron stars (BNS) in circular 
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orbits and show that the use of these techniques results in a notable improvement of the 
overall quality of the simulation. HR not only suppresses the Hamiltonian constraint 
violation but also contributes to a more stable behavior in the total angular momentum 
of the system. MR contribution is mostly confined to quenching the momentum 
constraint violation. The BNS simulation runs for about two orbits before stopping. 
One of the main reasons for the degradation of the simulation quality is the use of 
a shift vector frozen to its initial value (/3-freeze). This choice, while very convenient 
when testing new algorithms, becomes inadequate as soon as the stars move appreciably 
from their initial position. Given that our simulations are performed in the frame that 
rotates with the binary, this occurs rather late in the run. The second cause is related 
to inappropriate boundary conditions for the rest of the gravitational fields. The use 
of radiation (Sommerfeld) conditions in rotating frames becomes troublesome for large 
grids in rotating frames. These problems will be addressed in future work. 

Section El describes in detail the momentum constraint relaxation method and its 
numerical implementation. Section El presents simulations of BNS in circular orbits and 
compares the results obtained with and without the relaxation techniques. Section 0] 
summarizes our results and |Appendix A presents the convergence tests. 



2. Equations for the Gravitational and Hydrodynamical Fields 

2.1. Time Evolution Equations 

We use geometrized units (G — c — 1) and the Greek (Latin) indices run from to 3 (1 
to 3). In the standard "3+1" form, the metric is written as 

ds 2 = -a 2 dt 2 + j tj (da? + (3 i dt){dx j + ftdt) , 

where a, f3\ and ^ are the lapse function, shift vector, and spatial metric tensor, 
respectively. 

The ADM formulation [9 splits the Einstein's field equations 

into a set of time-independent elliptic equations 

R - KijK ij + K 2 = lQnp , (1) 
DjKj - DiK = 8nSi , (2) 

known as the Hamiltonian and momentum constraints, plus a set of time-dependent 
hyperbolic equations: 

(d t - Cp)jij = -2aK i:j , (3) 

(d t - CplKy = -DiDja + a{R tj - 2K il K l j + KK {j - 8tt[^ + -^{p - S)}} . (4) 

The latter set provides the evolution in time of the spatial metric 7^ and the extrinsic 
curvature K^. The symbol Di represents the covariant gradient with respect to the 
tensor 7^. The fields p, S, and are derived from the matter fields by splitting the 
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stress-energy tensor T^ v in components parallel and perpendicular to the normal of the 
spatial hypersurface n a |18,. 

Following York [EH EDI, we can decompose the tensors 7^ and Kij as 

7ij Tij : 

Where the fields if], 7^, A^-, and X, are known as the conformal factor, the conformal 
metric, the conformal traceless extrinsic curvature, and the trace of the extrinsic 
curvature respectively. We can write the Hamiltonian and momentum constraints using 
the new variables as 

fWiDjip -\r+ y A v Aij - %i K2 + 2tc ^p = °> ( 5 ) 

& (^Aji) - ^DiK - Sir^Si = . (6) 

We define the Hamiltonian constraint residual 7i and momentum constraint residual 
M.i as the l.h.s. of equations © and (jOJ) respectively. 

The BSSN formulation (TUJ HI] provides a set of evolution equations for the fields 

'ip, K, and 

(d t -C f3 )\n(^) = ~aK , (7) 

(d t - Cp)^ = -2ai 4j , (8) 

(d t - £p)K = -''DjDm + \aK 2 + aA^A^ + Ana{p + S) , (9) 

{dt - Cfs)A tJ = ^[-DiDja + a(R tJ - 8nS t:i )} TF + a{KA tj - 2A il A l j ) , (10) 

where the superscript TF indicates the trace-free part of the tensor. These fields are 
complemented with the variable known as the conformal connection 

r = -jf d , 

and its corresponding evolution equation 

= f k P\jk + \l ij P k ,ki ~ V >"., + + '"I"., - 2A i 1d j a 

- 2a (^Kj - 6i« KVOL" - ^M Jk + ^l ij Sj) . (11) 

ADM and BSSN are, in their original forms, unconstrained formulations. The 
Hamiltonian and momentum constraints are not enforced throughout the simulation but 
only monitored for quality control. Anderson & Matzner presented a constrained 
variation of the ADM formalism where the constraints are solved numerically at every 
time step. They use equation (J3J) to solve numerically for the conformal factor ip, while 
the momentum constraint equations (jBJ) are satisfied by the addition of a longitudinal 
component to the traceless conformal extrinsic curvature A l i as originally outlined by 
York - W e introduced in Paper I the Hamiltonian relaxation technique that draws 
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on the same idea of controlling the Hamiltonian constraint violating modes by adjusting 
ip. HR uses the conformal factor to solve approximately the Hamiltonian constraint (0), 
bypassing altogether the corresponding evolution equation (J7J for %j) . Instead of solving 
equation (J3J), HR composes an alternative parabolic equation which relaxes ip towards a 
solution of the Hamiltonian constraint through an iterative scheme. The equation used 
in Paper I was of the form 

d t tP = e H (d t H + VH H) , (12) 

where en and tjh are fine-tuning parameters. This technique efficiently quenches the 
development of Hamiltonian constraint violation instabilities, improving the overall 
quality of BNS simulations. 

In this paper we present the momentum constraint relaxation method, which 
suppresses the development of momentum constraint violation instabilities and works 
as a complement to HR. MR is based on correcting the traceless conformal extrinsic 
curvature A\j with a longitudinal component (Iw)^ [20]. The operator (Iw)^ is defined 
as 

2 

(Iw)^ = DiWj + DjWi - - jij D n w n , (13) 

and Wi is known as the vector potential. The momentum constraint equations 
become now a function of Wi 

M t = f 3 [ D k (lw) l3 + 6 (lw) l3 D k lnV] +Pi — 0, (14) 

where collects all terms independent of the vector potential 

D k Aij + 6 Aij D k In ^] - ^ - 8 tt S t . (15) 

Equations (|T4*jl form a set of coupled elliptic equations for the components of the 
vector potential w% that can be solved numerically using a number of different algorithms. 
Note a formal difference between HR and MR: HR replaces the BSSN evolution equation 
for the conformal factor (J7J) with a new equation (|T2*j) , while MR adds three new auxiliary 
fields, the Wi components, with their corresponding equations ((HJ). 

The components of the vector potential Wi and the conformal factor -0 are added 
to the list of dynamical fields evolved using BSSN: 7^, A^, K, and P. It is important 
to note that the fields Ay are not updated using the Wi new values. Instead, the uii are 
kept as dynamical variables that evolve from one time step to the next. The r.h.s. of the 
BSSN evolution equations are modified to include the contribution of u>i] in equations 
(jBHUJ), A^ is replaced by + This facilitates the implementation of Wi boundary 

conditions, as explained in section 12.3. 11 

The work presented in this paper is based on a Successive Under Relaxation 
algorithm (SUR) similar to the one described in [27]. This algorithm updates the value 
of Wi at each relaxation step r with a correction term 

— wl — uj m Awl > (16) 

where % is a relaxation parameter such that < uj m < 2 [27 \. When % > 1 
(i.e., Successive Over Relaxation or SOR) the method converges faster to the numerical 



Pi = f 3 
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solution than for the case lom < 1. However, a key concept behind momentum relaxation 
is to use the vector potential Wi to gradually project the simulation onto the space of 
solutions of the momentum constraint equations. Since we are interested in gentle 
updates of it>j, we use values of uj m < 1. SUR also allows for a rather straightforward 
implementation of the boundary conditions for the vector potential described in section 

mm 

Their values are updated at each one of the steps of the Iterative Crank-Nicholson 
(ICN) method (one Predictor and two Corrector steps) used for time integration. The 
implementation of momentum relaxation, which in this paper is used in combination 
with HR, is as follows : 

1) Initial update of Wi : After the initial data set corresponding to a BNS system 
in circular orbit is read, Wi undergoes a relatively large number of relaxation steps 
(typically about 90), starting from the initial guess w l = 0. The components are relaxed 
sequentially in the order w z , w y , and w x . Given that the three equations are coupled 
in the components of Wi, the last component to be relaxed experiences the smoother 
update. We choose the x component last because the violation of this component of the 
momentum constraint is the largest due to our choice of coordinate axis (binary aligned 
along the y axis and z = the orbital plane). 

2) At the Predictor stage of the ICN integration, we perform the update of the 
fields in the following order: 

i) Update of the gravitational fields 7^, Aij, K, and P. 

ii) Update of the conformal factor ip using HR (Paper I). 

iii) Update of the vector potential components w z , w y , and w x . The sources of 
the elliptic equations are calculated using the values of the fields corresponding to the 
previous time step. Each equation undergoes typically from 10 to 20 relaxation steps. 
The boundary conditions for Wi are updated right after each relaxation step (section 

iv) Update of the hydrodynamical fields, the lapse function, and the shift vector. 

The boundary conditions are updated at the same time as their corresponding 
fields. The same steps are followed in the first Corrector (second Corrector) stage, but 
replacing the field values from the previous time step with the corresponding updates 
generated in the Predictor (first Corrector) step. 

2.2. Lapse and Shift Equations 

We use the same lapse function and shift vectors employed in the development of HR 
(Paper I): the K-Driver algorithm j2E] for a and the shift vector remains unchanged 
from its initial value (/3-freezing). All the simulations presented in this paper were 
performed in the frame that rotates with the binary (corotating), given its superior 
stability properties over inertial frames 0injl2E|- We refer the reader to Paper I for the 
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details of the numerical implementation as well as the boundary conditions. 



2.3. Boundary Conditions 

We adopt Sommerfeld (radiative) boundary conditions for the conformal metric 7^, the 
traceless part of the extrinsic curvature Aij, and the trace of the extrinsic curvature K, 
and we set P = at the boundaries The boundary conditions for the conformal 
factor ip are such that enforce the satisfaction of TC = in its finite-difference form. 
They are essential for the control of the Hamiltonian constraint violating modes coming 
from the grid boundaries. 



2.3.1. Boundary Conditions for the Vector Potential W{ The boundary conditions for 
the vector potential Wi are inspired in their conformal factor counterparts and also 
aim at controlling as effectively as possible the incoming constraint violating modes. 
The boundary conditions for the vector potential iUj are such that they enforce the 
satisfaction to round-off error of the momentum constraint at the grid points next to 
the boundaries. Equation (jT4*)l can be expanded into 

Dk(lw)ij + 6 (lw) i3 , D k lnif) + Pi = 



D. D.u-, - D.Djir, - £ A D n w n + 6 f j 4(ln</>) 



2_ ~ 

DiWj + DjWi - -7y D n w r , 



j 1 ^h^j' 

+ Pi = 0, (17) 

where the first and second covariant derivatives of Wi are, as usual, functions of the 
partial derivatives with respect to the spatial coordinates. These partial derivatives of 
the vector potential are approximated in our finite-difference scheme by second order 
stencils of the form 



d x w? ~ 



d x d x w[ 



d x d y Wi ~ 
etc, 



w i(N,j,k) W i(N-2,j,k) 

W i(N,j,k) ~ ^ W ?(N-l,j,k) + W ?(N-2,j,k) 

(Kxf ' 

w i(N,j+l,k) ~ w i(N-2,j+l,k) ~ W i(N,j-l,k) + W i(N-2,j-l,k) 

4 Ax Ay 



where n is the relaxation iteration step, Ax and Ay the grid spacing along the x and y 
axis, respectively. The above formulae have been written for the grid points next to the 
boundary x = x max (i.e.; the points with indices (N — l,j,k)). When we replace the 
differential operators of equation (J 17)) with the finite-difference formulae (|18|) and solve 
for the boundary values of Wi at the boundary x = x max , we obtain 

W i(N,j,k) = F i(N,j,k) , (19) 

where Fi^,j,k) are algebraic expressions which depend on w^ N _ ± j^, vj^/ N j_ 1)f A, 
w i(N,j+i,k)i w i(N,j,k-i)i anc ^ w i(N,j,k+i)- The linear system of equations (0 is quite large 
even for modest-size grids and thus very time consuming to solve for. To expedite 
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the process, we evaluate the equations (|T^|) using the newly calculated values at the 

boundaries as soon as they are ready (i.e., "W-^tv-ij,*:)' w i(N,j-i,ky anc ^ w i(N,j,k-i)) an< ^ 
use the previous iteration values for the remainder (i.e., wf^j +1 ^ and ■ fe+1 ))- This 
approximation is justified by the fact that these boundaries are evaluated many times 
in every time step. 

Similar to the case of the boundary conditions for the conformal factor, a problem 
arises at the edges of the grid, where the expression for the boundary value wf, N N fe > is 
now a function of the inverse of the spatial conformal metric j xy , which is null at the 
initial time step. For these points we use bilinear interpolation until the value of % y 
crosses some predetermined value which, for the simulations of this paper, is 10~ 6 . 

2.4- Hydrodynamical Formulation and Initial Data sets 

The numerical simulation of neutron star spacetimes requires algorithms for the 
evolution of the matter fields. The fluid in the stars is described by a perfect-fluid 
stress-energy tensor and a polytropic equation of state with constant T = 2 is assumed. 
We use the hydrodynamical methods described in Paper I. The initial state of the BNS in 
circular orbit is obtained using the Wilson-Mathews Conformally Flat Condition (CFC) 
approach [2111201, using the elliptic solver described in [3T| . 

We tested the numerical implementation of the momentum relaxation method using 
short trial-and-error runs that consisted essentially in BNS simulations performed in 
small low-resolution grids. To simplify even further these runs, they were based on an 
initial data set corresponding to a corotating (tidally-locked) binary. The satisfactory 
performance of these short runs was followed by longer runs based on larger grids and 
higher resolution. These long runs were also based on more realistic irrotational (zero- 
spin) BNS systems. The details of the initial data sets used in these runs is given in 
table I of Paper I. The number of points in the short run grid is about 73 times smaller 
than that of the long run simulations, making it possible to simulate a BNS orbital 
period in a couple of hours on a typical single-processor workstation. The long runs 
were performed using the IA-64 Linux Beowulf cluster Mercury at NCSA. 

Our simulations are performed in Cartesian grids and use finite-difference second 
order operators within grids that have uniform and identical spacing along each axis. 
We work in the reference frame that rotates with the binary, and the stars are aligned 
with the y axis. 

3. Results 

3.1. Short Trial- and- Error Runs 

The free parameters of MR are the relaxation parameter lom and the number of 
relaxation steps. Their corresponding values were determined empirically and set to 
= 0.1 with 20 relaxation steps for the short runs. Note that ip and Wi are relaxed 
three times at every time step (i.e., in the ICN predictor and the two corrector stages). 
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Figure 1. Short Runs : Evolution of the three components of the momentum 
constraint residual during the first steps of the time evolution. The upper (lower) 
plot corresponds to a run without (with) momentum relaxation. The thick line is the 
initial residual. The curves are plotted following the line with coordinates (0, y,0), 
that runs through the center of the star. The companion star, located on the negative 
y hemisphere, is not shown. 
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Following the HR method used for the short runs of Paper I, we use the parameter 
values €h = 0.0001 and r/ H = 70.0, and the relaxation of if) is performed until the L 2 
norm of the Hamiltonian residual is smaller than the one at the previous time step for 
up to a maximum of 25 iterations. 

For comparison, we show in this and the next section curves corresponding to runs 
performed using only the BSSN formulation as described in [2], runs using only HR, and 
runs using both HR and MR. The BSSN runs differ from the HR and HR+MR runs in 
two aspects. One is that the BSSN runs use the T-Driver pffi lH] while the HR only and 
HR+MR runs use the /3-freeze condition for the shift vector. The other difference is that 
we use fall-off Robin-like boundary conditions for the lapse function for the BSSN runs 
and frozen values for the HR and HR+MR runs (see Paper I). All the runs (short and 
long) have Courant factors of 0.46. We also conserve the K-Driver and T-Driver (BSSN 
runs only) parameters, which are set to e a = 0.125, r/ a = 0.1, = 0.0005, r\p = 0.2, 
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Figure 2. Short Runs : Total angular momentum J as a function of time, given as 
fraction of the orbital period P. J is normalized to its initial value Jo- The lines 
correspond to runs using BSSN (dotted), HR only (thin solid), and HR + MR (thick 
solid). 
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respectively. The K-Driver (T-Driver) relaxation is iterated 5 (10) times. 

Figure presents the evolution of the momentum constraint residuals Aii for the 
first three time steps of a BNS simulation without (with) MR in the upper (lower) part 
of the plots. The momentum constraint violation is plotted along the binary axis (y). 
Only one star is present in the grid (the companion star is in the y < hemisphere) and 
its center is at y ~ 1.4. The thick lines correspond in all cases to the initial data. The 
reduction of the momentum constraint violation caused by MR can be seen both in the 
bulk of the grid and at the boundary y = y max . Note, however, that this effect is not 
as dramatic as in the case of HR and the Hamiltonian constraint violation reduction 
(compare with figures 1 and 2 of Paper I). Figure El shows the evolution of the total 
angular momentum of the BNS normalized to its initial value. Note that the use of MR 
(thick solid line) improves the quality of the run only marginally over the performance 
obtained with HR alone (thin solid line). The small grid size for which and orbital 
period is roughly equivalent to 10 side-to-side light crossing times makes this result 
quite remarkable. The BSSN simulation is presented with a dotted line. 
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Figure 3. Long Runs: Evolution of the Li norm of the momentum constraint 
violation across the numerical grid. The solid line corresponds to the MR+HR runs, 
while the dashed curve shows the result of using only HR. The improvement on the 
satisfaction of the momentum constraint is about a factor of 4 at the end of the 
simulation. 




3.2. Long High- Resolution Large-Grid Runs 

We tested the MR scheme with a large-grid high-resolution simulation of an irrotational 
BNS system. The details of the initial data set are provided in table 1 of Paper I and 
the corresponding convergence tests are presented in | Appendix A| All the plots of this 
section show curves that, for clarity, have been normalized to their corresponding initial 
values. The MR free parameter values used for these long runs were ujm = 0.01 with 20 
relaxation steps. 

The evolution of the L<i norm of the three components of the momentum constraint 
residual M.i across the numerical grid is shown in figure H3 While MR achieves the 
reduction of the constraint violation as it was designed to do, the improvement is not 
as impressive as in the case of HR and the Hamiltonian constraint violation (figure 4 of 
Paper I). At the end of the simulation, the momentum constraint violation was about 
four times smaller than in the HR only runs. The spikes present in the HR only curves at 
t ~ 0.4P occur on the stellar surface and are related to matter displacement in the grid, 
a side-effect of using a frozen shift vector. Note, however, that those spikes disappear 
when using MR. 

Figure |U shows the evolution of the total angular momentum. The plot shows the 
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Figure 4. Long Runs: Evolution of the total angular momentum J. The solid, 
dashed, dash-dotted lines corresponds to the MR+HR, HR only, and BSSN runs 
respectively. The dotted line shows the PN estimation (see Appendix B of Paper 
I). The inset expands the plot for the first half orbital period. 



1.15 



1.10 



1.05 



1.00 



0.95- 



0.0 









— \ p 




I : — 

/ 












l J 






i i 1 i i 


















rj 


_ 1.00 










1 
























1 


























- 0.99 
















i 1 i 1 i 1 i 1 i 






i 





0.1 0.2 0.3 0.4 0.5 




i 
































































i_ 











0.5 



1.0 
t/P 



1.5 



2.0 



MR+HR (solid), HR only (dashed) and BSSN (dashed-dotted) curves, as well as the 
PN estimation (dotted line) of the angular momentum loss for a point-mass binary with 
the same mass and angular momentum as the BNS in consideration (see Appendix B of 
Paper I). Both the HR only and MR+HR runs agree with the PN prediction for about 
1.5 orbital periods. The inset of figure 0] zooms in on the first half of the period, showing 
the reduced level of noise in the HR only and MR+HR curves. 

The contour plots |S] show in more detail the evolution of the momentum constraint 
violation during the simulation. The left (right) column shows snapshots of the BSSN 
(MR+HR) run at t=0, 0.5P, 1.0P, and 1.5P. The lower surface plots show the rest mass 
density, highlighting the position of the star in the grid. The instability that eventually 
stops this simulation originates at the corner of the cubical grid and is related to the 
use of radiation boundary conditions in combination with a rotating frame. While MR 
manages to reduce the effect of this constraint violation, simulations based on larger 
grids will have to address this problem. This effect is less serious in smaller grid runs, 
where the effect of frame rotation at the corner is consequently smaller. The analysis 
of this instability as well as a more comprehensive study of shift vectors for use with 
constraint relaxation methods will be done in future work. 

Figure El shows the remaining quality control curves: the coordinate separation 
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Figure 5. Long Runs: Surface plot of the violation of the x component of 
the momentum constraint. The left (right) column shows snapshots of the BSSN 
(MR+HR) run at t=0.5P, 1.0P, and 1.5P. The lower surface plots show the rest mass 
density, indicating the position of the star in the grid. 
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Figure 6. Long Runs: Remaining quality control curves for the MR+HR (solid), 
HR only (dashed), and BSSN (dash-dotted) runs. From top to bottom, we show the 
evolution of the coordinate separation between stellar centers d, the total gravitational 
M mass, and the L 2 of the Hamiltonian constraint. 
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between stellar centers d, the total gravitational M mass, and the L 2 of the Hamiltonian 
constraint. The latter curve corresponding to the BSSN run is not shown, given that it 
is out of scale. The total rest mass of the system remains invariant to within a 0.1 % in 
all runs. 

4. Conclusions 

We introduced the momentum relaxation method which complements the Hamiltonian 
relaxation scheme presented in Paper I. These algorithms gently relax the conformal 
factor if; and the vector potential Wi to achieve a gradual control on the growth of 
constraint violating modes. This smooth updating of the ip and Wi is a key ingredient 
of the methods since fully enforcing the satisfaction of the constraints (i.e. by relaxing 
the fields to obtain numerical solutions of the constraint equations) quickly renders the 
simulation unstable. 

The constraint relaxation methods have been tested in combination with BSSN 
in the simulation of binary neutron stars. Their use not only effectively quenches the 
constraint violating modes, but also improves the overall quality of the simulation as 
seen in the behavior of the total gravitational mass and angular momentum of the 
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binary. The simulations in this paper end after a couple of orbits due to two principal 
reasons: the inadequate use of the frozen shift condition and instabilities generated by 
the use of radiation (Sommerfeld) boundary conditions for the rest of the gravitational 
fields. The boundary conditions produce constraint violating modes at the corners of the 
cubical grid which grow worse with increasing grid size. Future work will concentrate on 
finding the best gauge choices as well as a way to avoid the problems generated at the 
corners of the grid. Another important aspect that remains to be studied is the behavior 
of the relaxation techniques in highly dynamical spacetimes. This will be addressed by 
studying their use in the simulation of BNS mergers. Finally, it remains to be seeing how 
these methods will affect the simulation in the presence of black holes. One potential 
complication could arise when using excision, since this would require the development 
of inner boundary conditions for the conformal factor and the vector potential. Note, 
however, that new techniques for black hole evolutions without excision have recently 
been developed [33J HHj • These methods have been succesfully employed in numerical 
simulations of black hole binaries [SHI OH E3 EE] and may prove to be a more robust 
platform for the implementation of relaxation techniques. 
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Appendix A. Code Tests 

To test the convergence of the results obtained when using the HR and MR, we 
performed the long run simulation on three different grids with lengths (in units of 
total rest mass) = 11.6, 14.0, and 18.6, while keeping the same grid spatial resolution 
( 40 grid point across the star). Figure IXT1 shows the evolution of the total angular 
momentum as a function of time. We see the convergence of the numerical results 
towards the PN estimation for point-mass binaries (dotted line). We also notice that 
the corner instability occurs sooner in the largest grid (solid line); the small grid run 
(dash-dotted line) stops due to the inadequacy of /3-freeze as a shift vector condition in 
the presence of matter displacement. 

To test the numerical convergence of the relaxation results with the spatial grid 
resolution, we performed three runs based on the irrotational long run using 20 (low 
res.), 25 (medium res.), and 40 (high res.) grid points across the star. However, in order 
to see explicitly the second order convergence we fully relaxed the vector potential Wi to 
achieve numerical solutions of the momentum constraint equations (jHJ). The results are 
shown in figure lA"2l One important effect is that the simulations become quickly unstable 
when using full relaxation in combination with BSSN, the K-Driver lapse function and 
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Figure Al. Convergence of the momentum relaxation results with varying grid 
sizes. The convergence test is based on the long run simulation and the plot shows 
the behavior of the total angular momentum. The labels for curves are (from smallest 
grid to largest) dash-dotted, dashed, and solid. All the grids have the same spatial 
resolution. The dotted line shows the PN estimation. 
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frozen shift vector. In all cases, the runs do not last beyond one tenth of an orbit. The 
reasons for this incompatibility between BSSN and the full enforcement of constraint 
satisfaction at every time step are not clear and deserve further study. 

The Hamiltonian relaxation technique uses an equation derived from the 
Hamiltonian constraint (J12J) to replace the BSSN evolution equation for the conformal 
factor (jZJ). Since the latter equation is not enforced anymore, we tested its satisfaction 
by evaluating the L 2 norm of its violation. We define the violation as 

\\d t ^BssN\\2 = ||<9tln(<0) - Cpki{ifi) + \aK\\ 2 

o 

(A.l) 

We implemented the finite-difference version of this equation approximating the time 
derivative with a second order central difference stencil centered at the time step t — 1 

HVbssnW* - \\ ln ^Y-^Y~ 2 _ (CM^r 1 + \{<*KT% . (A.2) 

Similarly, we tested the violation of the xx component of the traceless extrinsic curvature 
after the vector potential W{ is updated. Calculating A^ w as A^ w = A xx + (lw) xx , we 
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Figure A2. Evolution of the L2 norm of the violation of the components of the 
momentum constraint for the irrotational Long Run described in table I of Paper I 
with three different grid resolutions. The dotted, dashed, and solid curves correspond 
to the low, medium, and high resolution runs and the numerical factors multiplying 
the curves correspond to the ratios between grid spacings. 
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define the violation as 

ANew t _ ANew (t-2) 

\\dtA»: w BSSN \\ 2 « ||^ ^Af i^ry- 1 - {[rhs A^}y-% , (A.3) 

where [rhs A^ w ] represents the r.h.s. of equation (|10|) . 

Figure IA3I shows the behavior of such violations for the small runs described in 
section EHJ The comparison of H^t^BSSJvlh (top) from the HR+MR run (solid curve) 
with the corresponding violation from the BSSN run (dashed curve) shows no significant 
difference |. The plot of \ \dtA^ w BSSN \\2 (bottom) reflects a drift in the BSSN curve that 
is not present in the HR+MR case. This apparent "improvement" due to the use of 
relaxation methods is likely due to the overall stability gained throughout the simulation 
and not a specific betterment of the satisfaction of equation (jlOj) . 



| Note that the BSSN curve is not identically zero since the second order convergence in time is achieved 
in the evolution through the use of ICN integration and not the time derivative stencil of eauation lA.2l 
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Figure A3. Evolution of the L2 norm of the violation of the BSSN evolution equation 
for the conformal factor (top) and the xx component of the traceless extrinsic curvature 
(bottom). The dashed and solid curves correspond to the BSSN and HR+MR short 
runs described in section l3~Tl 
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